Transcriptome-wide map of N6-methyladenosine (m6A) profiling in coronary artery disease (CAD) with clopidogrel resistance

Background Clopidogrel resistance profoundly increases the risk of major cardiovascular events in coronary artery disease (CAD) patients. Here, we comprehensively analyse global m6A modification alterations in clopidogrel-resistant (CR) and non-CR patients. Methods After RNA isolation, the RNA transcriptome expression (lncRNA, circRNA, and mRNA) was analysed via RNA-seq, and m6A peaks were identified by MeRIP-seq. The altered m6A methylation sites on mRNAs, lncRNAs, and circRNAs were identified, and then, GO and KEGG pathway analyses were performed. Through joint analysis with RNA-seq and MeRIP-seq data, differentially expressed mRNAs harbouring differentially methylated sites were identified. The changes in m6A regulator levels and the abundance of differentially methylated sites were measured by RT-PCR. The identification of m6A-modified RNAs was verified by m6A-IP-qPCR. Results The expression of 2919 hypermethylated and 2519 hypomethylated mRNAs, 192 hypermethylated and 391 hypomethylated lncRNAs, and 375 hypermethylated and 546 hypomethylated circRNAs was shown to be altered in CR patients. The m6A peaks related to CR indicated lower mark density at the CDS region. Functional enrichment analysis revealed that inflammatory pathways and insulin signalling pathways might be involved in the pathological processes underlying CR. The expression of mRNAs (ST5, KDM6B, GLB1L2, and LSM14B), lncRNAs (MSTRG.13776.1 and ENST00000627981.1), and circRNAs (hsa_circ_0070675_CBC1, hsa-circRNA13011-5_CBC1, and hsa-circRNA6406-3_CBC1) was upregulated in CR patients, while the expression of mRNAs (RPS16 and CREG1), lncRNAs (MSTRG.9215.1), and circRNAs (hsa_circ_0082972_CBC1) was downregulated in CR patients. Moreover, m6A regulators (FTO, YTHDF3, and WTAP) were also differentially expressed. An additional combined analysis of gene expression and m6A peaks revealed that the expression of mRNAs (such as ST5, LYPD2, and RPS16 mRNAs) was significantly altered in the CR patients. Conclusion The expression of m6A regulators, the RNA transcriptome, and the m6A landscape was altered in CR patients. These findings reveal epitranscriptomic regulation in CR patients, which might be novel therapeutic targets in future. Supplementary Information The online version contains supplementary material available at 10.1186/s13148-023-01602-w.


Introduction
Clopidogrel, a novel antiplatelet drug used to treat CAD patients after percutaneous coronary intervention (PCI), reduces platelet aggregation by inhibiting the activity of the adenosine diphosphate (ADP) receptor P2Y 12 .However, approximately 10-30% of treated patients still suffer from ischaemic events [1], which might be due to a reduced degree of platelet inhibition.The inability to effectively prevent platelet aggregation, which results in high residual platelet activity, is considered to be indicative of low response to clopidogrel or clopidogrel resistance.It clearly increases the risk of recurrent cardiovascular and cerebrovascular events in CAD patients after PCI.
Complex factors have been reported to influence the clopidogrel response, including clinical factors (such as diabetes), sequence alterations (SNPs, lncRNAs, miR-NAs, and mRNAs), and epigenetic modifications (DNA methylation and RNA methylation).Among these factors, RNA methylation is currently a hot research topic.In most eukaryotes, m6A modification in the 5′ cap region of mRNA plays an important role in the maintenance of mRNA stability, variant splicing, polyadenylation, transportation, and translation, and a m6A mark in the 3′ polyA region contributes to the extracellular transport, initial translation, and structural stability of the mRNA along with polyA-binding proteins [2].The dynamic modification of m6A methylation is related to life processes such as immune responses and stem cell renewal in the body, which can affect the signalling mediated through downstream molecules, cause gene expression imbalance, and alter cell differentiation, homeostasis, and the stress response, leading to the occurrence of diseases or acquisition of pathological states.
However, few studies have reported the effects of m6A modification on the efficacy of cardiovascular drugs.The previous studies have confirmed that METTL3-mediated m6A modification upregulated the expression of the long-chain intergenic noncoding RNA linc00958, leading to increased adipogenesis and affecting the prognosis of hepatocellular carcinoma chemotherapy [3].The total saponin content in Panax notoginseng can affect the WTAP/p16 signalling pathway through FTO-mediated m6A modification, thereby inhibiting the proliferation, migration, and intimal hyperplasia of vascular smooth muscle cells [4].However, whether this dynamic modification can affect the antiplatelet effect of clopidogrel is unknown.In the present study, our group sequenced the whole transcriptome from CAD patients, and then, the m6A methylation abundance of the whole transcriptome was analysed.An integrative analysis of RNAs revealed the molecular mechanism underlying clopidogrel resistance.

Characteristics of the patients
There were 46 patients enrolled in our study, and samples from ten patients (five CR and five NC) were used for the systematic analysis of RNA m6A modification abundance, and samples from the remainder of the patients (18 CR and 18 non-CR) were used for validation.The clinical features of the participants are shown in Table 1.Except for that of the platelet function measure, no significant difference in baseline levels was found between the two groups.
A total of 18 CR and 18 non-CR patients with CAD were analysed for the correlation study.The baseline values of the parameters for these patients relisted in Table 2.In addition to the platelet function index and BMI values, the characteristics of the patients in the validation cohort were closely matched.

RNA transcriptome expression profiles in CR and non-CR patients
To further explore the relationship between gene expression and m6A modification, we measured the levels of potential markers of heterogeneity after clopidogrel treatment.We characterized the transcriptome in CR and non-CR patients by RNA-seq with high throughput on the basis of circRNAs, lncRNAs, and mRNAs.With thresholds set to |log2FC|> 1 and a P value < 0.05, we identified 583 lncRNAs (391 downregulated and 192 upregulated) and 922 mRNAs (546 downregulated and 376 upregulated).In the light of the expression profile of the circRNA microarray, long probes and short probes were used, and 5448 dysregulated circRNAs (2529 downregulated and 2919 upregulated) were identified at threshold levels of an |log 2 FC|> 1 and a P value < 0.05 (Additional file 1: Stable 1, Additional file 2: Stable 2, and Additional file 3: Stable 3).

Overall characteristics of the m6A modification in CR and non-CR patients
Blood samples from ten patients (five CR and five non-CR) were collected and used for MeRIP-seq analysis.As the Venn diagram in Fig. 1A-C displays, 38,071 m 6 A peaks on mRNAs, 21,454 m 6 A peaks of lncRNAs, and 15,944 m 6 A peaks on circRNAs were discovered between CR and non-CR patients.Among these peaks, 751 nonoverlapping m6A peaks were found on mRNAs (Fig. 1A), 446 nonoverlapping m6A peaks were found on lncRNAs (Fig. 1B), and 78 nonoverlapping m6A peaks were found on circRNAs (Fig. 1C) in CR patients.In addition, 114 nonoverlapping m 6 A peaks were found on mRNAs (Fig. 1A), 152 nonoverlapping m 6 A peaks were found on lncRNAs (Fig. 1B), and 45 nonoverlapping m 6 A peaks were found on circRNAs (Fig. 1C) in the non-CR patients.
Then, we measured the distribution of m 6 A marks in the whole transcriptome.The vast majority of m6Amodified RNAs were covered by fewer than 5 m 6 A peaks, and a small number were covered by more than 5 peaks (Fig. 1E-F).
We used HOMER motif software to analyse the m 6 A peaks.In contrast with that in the CR groups, the motif sequence in the mRNAs was GGACC in the non-CR groups (Fig. 2A), which was consistent with the previously obtained results.To explore the potential function of m6A, we measured the distribution of m6A peaks across the whole transcriptome.We found that the m6A peaks were mostly enriched at the start of the 3′UTR in both the CR and non-CR patients samples (Fig. 2B).m6A peaks in the CR patient samples showed a lower density in the CDS region than those in the non-CR patient samples.In addition, the distribution of m6A peaks in the CR and non-CR patient samples showed had the same tendency, increasing from the 5′UTR to the start of the 3′UTR and decreasing from the start of the 3′UTR to the end of the transcriptome.We determined that the proportion of m6A peaks was the largest in the 5′UTR and the smallest in the intron in the non-CR patient sample (Fig. 2C and D).
These hypermethylated and hypomethylated m6A sites in mRNA, lncRNA, and circRNA were then classified by transcript region.As shown in Fig. 5D, the fold change was highest for introns with both hypermethylated and hypomethylated m6A sites in mRNA.For the hypermethylated sites in lncRNAs, the fold change was the highest downstream.For the hypomethylated sites in lncRNAs, the fold change was the highest in the promoter (Fig. 5E).Moreover, the fold change was the highest in the promoter for hypermethylated sites in circRNA, and then, it was the highest in the intron for hypomethylated sites in circRNA (Fig. 5F).From all these results, we concluded that there were distinct differences in m6A marks both in terms of their distribution and in their abundance between CR and non-CR patients, indicating that different m6A modification patterns may lead to clopidogrel resistance.

Functional annotation of altered m6A-modified mRNAs, lncRNAs, and circRNAs between CR and non-CR patients
To characterize the potential function of differential m6A methylation patterns between CR and non-CR patients, mRNAs, lncRNAs, and circRNAs with distinct m6A abundance levels were subjected to Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Gene and Genomes pathway analysis.
GO and KEGG enrichment analyses with circRNAassociated genes that were differentially methylated were performed.For the biological process category, altered m6A-modified genes were significantly interrelated with the regulation of the RNA biosynthetic process (Fig. 6A  and D).For the cell composition category, hypomethylated circRNA-associated genes were mainly enriched in platelet alpha granule membrane, cytosol, and nuclear  6A), while hypermethylated circRNA-associated genes were significantly associated with the nucleoplasm, nuclear lumen, and nuclear body (Fig. 6D).For the molecular function category, downregulated circRNAassociated genes were involved in ubiquitin-like protein ligase activity, ubiquitin protein ligase activity, and protein binding (Fig. 6A), while upregulated circRNAassociated genes were mainly enriched in nucleic acid binding, transcription coactivator activity, and transcription cofactor activity (Fig. 6D).KEGG pathway analysis indicated that circRNA-associated genes were enriched in inflammatory pathways (such as the MAPK, Wnt, and Notch signalling pathways) (Fig. 7A and B).
In Fig. 6B, for the BP, CC, and MF categories, the hypomethylated lncRNA-associated genes were enriched mainly in lipid metabolic process, glycosphingolipid metabolic process, and negative regulation of protein catabolic process; lysosomal lumen, vacuolar lumen, and primary lysosome; hydrolyzing O-glycosyl compounds, protein phosphatase binding, and acting on glycosyl bonds, respectively.The hypermethylated lncRNA-associated genes, as displayed in Fig. 6E, were enriched mainly in UV protection, cAMP catabolic process, and cyclic nucleotide catabolic process (BP); intracellular part, intracellular, and cytoplasm (CC); enzyme binding, protein domain-specific binding, and cAMP binding (MF).KEGG pathway analysis revealed that differentially methylated lncRNA-associated genes were also significantly related to inflammatory pathways (such as Notch and Toll-like receptor) (Fig. 7C  and D).
In Fig. 6C, for the biological process category, mRNAs harboured m6A sites with downregulated methylation were involved in cotranslational protein targeting to the membrane, SRP-dependent cotranslational protein targeting to the membrane, and protein targeting to the ER.However, m6A sites with upregulated methylation were highly involved in nucleic acid metabolic processes, transcription, DNA templating, and RNA metabolic processes (Fig. 6F).For CC terms, downregulated m6A site methylation was significantly enriched in cytosolic ribosome, cytosolic part, and ribosomal subunit (Fig. 6C), while upregulated m6A site methylation was enriched in the nucleus, nuclear lumen, and nucleoplasm (Fig. 6F).Finally, these mRNAs were significantly correlated with nucleic acid binding, an MF category term (Fig. 6C and  F).Furthermore, these downregulated mRNAs were found to be significantly involved in the insulin signalling pathway, VEGF signalling pathway, and so on.The upregulated mRNAs were also involved in the insulin signalling pathway and some inflammatory pathways (such as Notch and mTOR) (Fig. 7E and F).

Comprehensive analysis of MeRIP-seq and mRNA-seq data in CR and non-CR patients
The 25 most highly differentially expressed mRNAs are displayed in Table 5.Through a comprehensive analysis of the data from RNA-seq and MeRIP-seq of CR and NCR patients, we identified one hypermethylated and upregulated (hyper-up) gene (LYPD2), one hypomethylated and downregulated (hypo-down) gene (RPS16), and two hypermethylated and downregulated (hyper-down) genes (KDM6B and ST5) with differential expression in CR patients compared to that in non-CR patients (Fig. 8).

Expression of m6A regulators and candidate genes correlated with clopidogrel resistance
In an effort to determine whether m6A methylation status was changed in CR patients, with the validated cohort data, we assessed m6A regulator expression in CR and non-CR patients through RT-PCR.The results indicated that the expression of FTO, a main m6A eraser (P = 0.0003), and YTHDF3, a main m6A reader (P = 0.0106), was significantly upregulated in CR patients (Fig. 9F and L) compared with non-CR patients (controls).However, the expression of WTAP, a main m6A writer, was lower (P = 0.0244) (Fig. 9D).The expression of other m6A regulators was unchanged between the CR and non-CR patient groups (Fig. 9A-C, E, and G-K).

Discussion
Antiplatelet treatment is the key treatment for patients with CAD after PCI.Although ticagrelor, a new antiplatelet drug, exhibits rapid, effective, and sustained inhibition of the P2Y12 receptor, patient haemorrhage risk is markedly increased [5].In clinical practice in China, although ticagrelor reduces cardiovascular mortality and the risk of stent thrombosis in ACS patients, it has an advantage only in patients with low bleeding risk and increases the incidence of cardiovascular events in patients with medium to high bleeding risk [6].The 2018 COSTIC study proposed that clopidogrel shows a significantly greater clinical benefit for the Chinese population [7].For patients with stable coronary heart disease or elderly patients, clopidogrel might be more beneficial.However, clopidogrel resistance significantly increases the risk of adverse cardiovascular and cerebrovascular events.The relevant mechanisms underlying the occurrence of clopidogrel resistance are complex.Various studies have indicated that the extent of clopidogrel resistance is influenced by a series of clinical morbidities, including chronic kidney disease, diabetes mellitus, and smoking and factors such as drug-drug interactions (for instance, proton-pump inhibitors), age, inflammatory status, decreased left ventricular ejection fraction, and increased body mass index [8].The molecular mechanism underlying clopidogrel resistance has been explored in multiple fields on the basis of sequence differences and epigenetic regulatory factors.First, based on singlenucleotide polymorphisms (SNPs) in drug metabolismrelated genes, researchers found that ABCB1 genes are involved in duodenal transport; CYP family genes are involved in liver metabolism; PON-1 genes and P2Y12 genes are related to platelet activation.Preliminary research by this project team revealed that the ABCB1 C3435T polymorphism is associated with platelet activity and MACEs in ACS patients taking clopidogrel [9].Another genome-wide association study suggested that genetic variations in genes such as CYP2C19, PEAR1, and N6AMT1 are associated with different antiplatelet effects and clinical endpoint events related to clopidogrel [10].Second, the up-or downregulation of noncoding RNA also affects the antiplatelet activity of clopidogrel.A recent omics study suggested that upregulation of the lncRNAs NONHSAT083775.2 and NONHSAT107804.2 or downregulation of the lncRNA NONHSAT133455.2 led to clopidogrel resistance [11].Another study with DM animal models revealed that lncRNA MT1P3 (metallothionein 1 pseudogene 3) promoted the expression of P2Y12 receptors by sponging miR-126, thereby increasing platelet activity [12].Moreover, miRNAs participate in the regulation of clopidogrel resistance [13]; for instance, miRNA-26a upregulation increased the clopidogrel resistance risk in patients undergoing PCI [14], and platelet-derived miRNAs (such as miRNA-24-3p, miRNA-142-3p, and miRNA-411-3p) might be potential biomarkers.In addition, the advancement of epigenetic studies has provided a new perspective on the clopidogrel resistance mechanism.Our group found that in CAD patients with hyperlipidaemia, CpG4 methylation of the PON1 promoter resulted in low mRNA expression of PON1, leading to clopidogrel resistance [15].Our latest research was based on a ChIP analysis of methylated DNA from clopidogrel resistance samples and revealed the identification of 979 hypermethylation sites and 6119 hypomethylation sites, and four sites (cg23371584, cg15971518, cg04481923, and cg22507406) were verified to be clopidogrel resistance diagnostic markers [16].In addition, N6-methyladenosine (m6A) might play a vital role in posttranscriptional modification.
In the present study, by preparing the RNA transcriptome expression profiles and validating the findings with samples from a larger cohort, we discovered that the expression of mRNAs (ST5, KDM6B, GLB1L2, and LSM14B), lncRNAs (MSTRG.13776.1 and ENST00000627981.1),and circRNAs (hsa_ circ_0070675_CBC1, hsa-circRNA13011-5_CBC1, and hsa-circRNA6406-3_CBC1) was upregulated in CR patients, while mRNAs (RPS16 and CREG1), lncRNAs (MSTRG.9215.1),and circRNAs (hsa_circ_0082972_ CBC1) were downregulated in CR patients.This is the first study to report on the landscape of the transcriptome based on sequencing of samples from patients with clopidogrel resistance, and it revealed some CRrelated biomarkers.Although the results need to be validated with data from a larger sample and from a broader population, the markers it revealed will provide a basis for early detection in clinical practice.Thus, the ST5, KDM6B, and CREG1 genes deserve our continued attention.First, the Cryptococcus neoformans sequence type 5 (ST5) lineage can be used to infect immunocompetent hosts and increase the blood platelet count [17].In our study, the ST5 level was higher in CR patients, which might have been due to a high level of ST5 increasing the platelet count, leading to clopidogrel resistance.Second, lysine-specific demethylase 6B (KDM6b) is involved in vascular remodelling, and it can promote neointima formation by inhibiting the activation of Nox4 signalling in autophagy [18].Third, the cellular repressor of E1A-stimulated genes (CREG1), which was downregulated in CR patients, promotes HUVEC proliferation and prevents endothelial cells (ECs) from undergoing apoptosis.Through ILK-Cdc42 activation, CREG1 can also increase EC filopodia formation and vascular assembly to promote neovascularization [19].Likewise, a latest study uncovered the role of CREG1 in the regulation of megakaryocyte maturation and thrombopoiesis [20].Hence, ST5, KDM6B, and CREG1 might be therapeutic targets for clopidogrel resistance and vascular diseases.Future researches might further validate the findings of this article.
Moreover, we carried out MeRIP-seq with samples from patients with low clopidogrel responses and found that m6A modification of mRNA on ST5, LYPD2, RPS16, PIP5K1C, RPS21, SRRT, and RPS6 affected mRNA expression, resulting in clopidogrel resistance.In addition, the m6A regulators FTO, YTHDF3, and WTAP were also differentially expressed in CR and non-CR patients.These findings may indicate that upregulated FTO or downregulated WTAP causes reduced expression of related genes.FTO, a key m6A enzyme, is involved in the pathophysiological processes leading to myocardial hypertrophy and cardiac dysfunction [21].Depending on the m6A modification, FTO regulates the cardiac function after systolic heart dysfunction in model rats and influences cardiac remodelling and repair [22].Loss of the m6A reader YTHDF2 promotes demethylation of H3K27me3, leading to enhanced production of proinflammatory cytokines [23].Moreover, WTAP has been reported to regulate the progression of myocardial ischaemia/reperfusion injury via the YTHDF1/FOXO3a signalling pathway [24].After our comprehensive analysis, we speculated that the higher degree of m6A methylation on ST5 mRNA might lower its expression and result in clopidogrel resistance.The intrinsic molecular mechanism needs to be confirmed with cell experiments.
Furthermore, in this study, we proposed that inflammatory pathways (such as the NF-κB, Wnt, and Notch signalling pathways) as well as the insulin signalling pathway participate in the regulation of clopidogrel resistance.For example, through the NF-κB signalling pathway, nucleotide-binding oligomerization domain 2 (NOD2) mediates P2Y12 upregulation and increases platelet activation, resulting in clopidogrel resistance [25].On the other hand, abnormal excretion of insulin leads to an increase in platelet activity, and hyperglycaemia activates the protein kinase C pathway by saccharifying proteins on the surface of platelets, inducing the expression of P selectin, reducing the fluidity of the cell membrane, and thus increasing the platelet adhesion rate and enhancing platelet activity [26].FTO is an important factor in blood glucose regulation, because it regulates insulin signallingrelated gene expression possibly through m6A modification.Research suggested that FTO is associated with insulin secretion and abnormal blood glucose levels, and it may interact with transcription factors (STAT3, CREB, ATF4, and FoxO1) and thus change their activity, regulating the expression of liver gluconeogenesis genes (PCK1 and G6PC) [27].A functional experiment with pancreatic ß cells showed that increasing the expression of FTO led to an increase in the gluconeogenesis rate and that reducing the activity of FTO led to the opposite results [28].Hence, FTO may regulate insulin signalling-related gene expression through m6A modification, changing the antiplatelet effect of clopidogrel.Further animal and cell experiments are needed to confirm this possibility.
To the best of our knowledge, our research is the first to report RNA transcriptome expression profiles and potential impact of m6A methylation on clopidogrel resistance.Although we made great efforts, this research still has unavoidable limitations.First, it was a single-centre study, and the studied population was not large.We will further expand the sample size to verify our conclusions.Second, PBMCs are collected from peripheral blood, and their status may be influenced by various factors, such as individual differences, etc., which may interfere with the interpretation of RNA transcriptome, and RNA methylation may be influenced by other factors, such as environment factors and coexisting diseases.These influences may have affected our results.Finally, knockout animal models need to be used for in vivo validation studies, and an in-depth mechanistic investigation will increase our understanding in the future.

Conclusions
In general, our present study revealed the RNA transcriptome expression profiles in CR and non-CR patients and showed that the expression of some mRNAs (ST5, KDM6B, GLB1L2, and LSM14B), lncRNAs (MSTRG.13776.1 and ENST00000627981.1),and circR-NAs (hsa_circ_0070675_CBC1, hsa-circRNA13011-5_ CBC1, and hsa-circRNA6406-3_CBC1) was upregulated in CR patients, while that of other mRNAs (RPS16 and CREG1), lncRNAs (MSTRG.9215.1),and circRNAs (hsa_circ_0082972_CBC1) was downregulated in CR patients.Moreover, m6A regulators (FTO, YTHDF3, and WTAP) were differentially expressed.By MeRIP-seq, we found that 2919 hypermethylated and 2519 hypomethylated mRNAs, 192 hypermethylated and 391 hypomethylated lncRNAs, and 375 hypermethylated and 546 hypomethylated circRNAs were markedly changed in CR patient samples.The joint analysis of mRNA m6A peaks and expression revealed that the mRNAs (such as ST5, LYPD2, and RPS16) might cause clopidogrel resistance because their expression is changed via m6A modification.These are critical regulators that interfere with the epigenetic regulation of clopidogrel resistance via the inflammatory pathway and insulin signalling pathway, and these findings provide new insights useful for CR mechanistic research and will help clinicians identify individualized treatments more effectively.

Patients
Between April and December 2020, a total of 23 CR and 23 non-CR controls in the First Affiliated Hospital of Ningbo University were enrolled in the present research.All the enrolled patients were diagnosed with ACS and were treated of PCI with drug-eluting stents to re-establish the coronary artery blood supply.The case group conformed to the standard for a CR diagnosis, which is reaction units of P2Y12 (PRU) greater than 240 [29] as measured by the VerifyNow P2Y12 assay (Accumetrics, Inc., San Diego, California).The enrolled CR and non-CR patients did not have a history of acute or severe infection, disordered kidney or hepatic function, active haemorrhage, or rheumatoid-related diseases.The clinical data, demographic data, and laboratory indices of these patients were collected.These data included BMI, age, sex, drinking rate, smoking rate, comorbidities (diabetes, hypertension, or hyperlipidaemia), medications, aTnI, ALT, AST, LDL-C, creatinine, uric acid, platelets, pro-BNP status, and other measures.

RNA preparation
Peripheral blood was collected from 23 CR patients and an equal number of non-CR subjects after clopidogrel therapy that was administered longer than 6 months.Samples from 10 of these patients were used for a comprehensive analysis of RNA m6A modification (the five pairs of CR and non-CR patients were well matched), while samples from the remainder of the patients (18 CR and 18 non-CR patients) were used to validate the results.
Peripheral samples were obtained from the patients described above and placed into an ETDA anticoagulant vacutainer.Peripheral blood mononuclear cells (PBMCs) were isolated as described previously [30].Total RNA was obtained with TRIzol reagent (Invitrogen, USA), and the RNA purity and concentration were determined with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA).

RT-PCR
RT-PCR was performed with cDNA.First, reverse transcription was performed with RNA by SuperScript III Reverse Transcriptase (Invitrogen, USA) according to the manufacturer's instructions.The following reagents were added in sequence to a 0.5-ml PCR tube: 2 µl of first-strand cDNA; 2 µl of upstream primer (10 pM); 2 µl of downstream primer (10 pM); 4 µl of dNTP (2 mM); 5 µl of 10 × PCR buffer; and 1 µl of Taq enzyme (2 u/ ul).Then, PCR was performed according to the following programme: 95 °C, 10 min; 30 PCR cycles (95 °C, 10 s and 60 °C, 60 s), and amplification reaction 28-32 cycles (95 °C, 10 s; 60 °C, 60 s; and 95 °C, 15 s) and slowly warmed from 60 to 99 °C to establish the melting curve.The housekeeping genes U6/GAPDH were utilized as the internal controls to ensure the reliability and accuracy of the experimental results.

RNA-seq and MeRIP-seq
RNA-seq and MeRIP-seq were performed at CloudSeq Biotech Inc. according to a procedure described previously [31].Briefly, after RNA isolation, the total RNA was fragmented into approximate 100 nt sequences, and those fragments were incubated with anti-m6A antibody (Manga) for 2 h at 4 °C.Then, total RNA was incubated with prepared beads for 2 h at 4 °C.Finally, the complexes were washed from the beads and purified with TE buffer.Both the RNA that was not immunoprecipitated and IP samples were used for library construction.The library quality was evaluated with a Bioanalyzer 2100 system.The library was sequenced on an Illumina HiSeq instrument.All the sequence data have been uploaded Supplement date (Additional file 4: Stable 4, Additional file 5: Stable 5, and Additional file 6: Stable 6).

Data analysis of MeRIP-seq and RNA-seq
Paired-end reads were obtained using a Illumina HiSeq 4000 sequencer, and the quality was controlled on the basis of a Q30 sore.Then, 3′ adaptor trimming and lowquality reads were removed by Cutadapt software (version: 1.9.3).Then, all the clean reads in the library were aligned to the reference genome (UCSC HG19) by Hisat2 software (version: 2.0.4).DCC software was used to recognize circRNAs based on the results of STAR alignment.Methylated sites on RNAs (m6A peaks) were recognized using MACs software.The recognized methylated sites were used for motif enrichment analysis by HOMER [32].In addition, the R package MetaPlotR was used to characterize the m6A distribution.diffReps was used to identify differentially methylated sites at threshold levels of |log 2 FC| higher than 1 and a P value less than 0.05.Scripts were written in-house to identify and select the differentially methylated sites that overlapped with mRNA exons.For RNA-seq, raw counts were obtained with HTSeq software (version 0.8.2), followed by normalization using edge R software.Genes of interest were directly visualized with Integrative Genomics Viewer (IGV; version: 2.3.68).

GO and KEGG pathway enrichment analysis
Gene Ontology (GO) enrichment analysis in three categories: biological process (BP), cellular component (CC), and molecular function (MF) categories, was performed to identify the function of identified genes related to differential peaks or differential expression of mRNAs, lncRNAs, and circRNAs (http:// www.geneo ntolo gy.org).GO analysis of differentially expressed peaks or RNAs was performed by using R software based on hypergeometric distribution.The number of genes related to altered peaks or enriched in each GO term was counted, and the related gene attributes annotated with a GO term were identified (http:// www.geneo ntolo gy.org).
In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (www.genome.jp/ kegg) was performed by using R software, and a hypergeometric distribution test was performed to evaluate the significance of genes related to altered peaks or that of differential expression of mRNAs, lncRNAs, and cir-cRNAs in association with each pathway term.

m6A-IP-qPCR and RT-qPCR
Twenty-two genes with differentially methylated sites as determined according to MeRIP-seq were verified by RT-qPCR.A small amount of fragmented RNA was used as the input control.The remainder of the RNA was incubated with anti-m6A antibody-coupled beads.The m6A-interacting RNAs were then immunoprecipitated and eluted from the beads.Both input control and m6A-IP samples were subjected to RT-qPCR with gene-specific primers.

Statistical analysis
All statistical analyses were performed using GraphPad Prism 8.0 software.The classification data were described as percentages, and the differences between individual groups were determined by Chi-square test.The measurement data were expressed as the mean ± standard deviation, and the significance of differences between two groups was determined by t-test or Pearson correlation test.A two-tailed P value less than 0.05 was considered to be statistically significant.

Fig. 1
Fig. 1 Overview of m6A modification peaks on mRNAs, lncRNAs, and circRNAs in CR and non-CR patients.A-C Venn diagram depicting the common and distinct m6A peaks in mRNAs, lncRNAs, and circRNAs between CR and non-CR patients.D-F The numbers of m6A peaks located on mRNA, lncRNA, and circRNA between CR and non-CR patients

Fig. 2 Fig. 3
Fig.2Overview of m6A peaks located on mRNAs, lncRNAs, and circRNAs in CR and non-CR patients.A Motifs enriched of all the identified m6A peaks in CR and non-CR patients.B Diversity of the m6A peak density in the indicated regions between CR and non-CR patients.C Pie chart revealing the regional distributions of m6A peaks in the RNA transcriptome of CR patients.D Pie chart revealing the regional distributions of m6A peaks in the RNA transcriptome of non-CR patients

Fig. 4
Fig. 4 Hierarchical clustering analysis showing the differences in the RNA transcriptome in CR and non-CR patients based on |log 2 FC|> 1 and P value < 0.05 criteria.A circRNAs, B lncRNAs, and C mRNAs.In heatmaps, red indicates hypermethylation, and green indicates hypomethylation

Fig. 5 Fig. 6
Fig. 5 Differences in the chromosomal distribution of differentially methylated m6A sites between CR and non-CR patients.A-C Chromosomal distribution of differentially methylated m6A sites in three different kinds of RNAs.D-F Statistical analysis showed fold change differences in hypermethylated and hypomethylated m6A sites in transcript regions of mRNA, lncRNA, and circRNA in CR versus non-CR patients

Fig. 7
Fig. 7 KEGG pathway analysis results showing differentially abundant m6A marks between CR and non-CR patients.A circRNAs in CR patients, B circRNAs in non-CR patients, C lncRNAs in CR patients, D lncRNAs in non-CR patients, E mRNAs in CR patients, and F mRNAs in non-CR patients

Fig. 8
Fig. 8 Comprehensive analysis of mRNA m6A peaks and their expression between CR and non-CR patients.Scatter plot showing the distribution of mRNAs with both significantly changed m6A and mRNA levels in CR and non-CR patients

Fig. 9
Fig. 9 The mRNA expression of m6A regulators in CR and non-CR patients.A-E m6A writers (METTL3, METTL14, METTL4, WTAP, and KIAA1429), F, G m6A erasers (FTO and ALKBH5), and H-L m6A readers (YTHDC1-2 as well as YTHDF1-3) This research was approved by the Ethics Committee of the First Affiliated Hospital of Ningbo University, and every individual signed a written informed consent form before the study was initiated.The study was performed in accordance with Declaration of Helsinki principles.

Fig. 10
Fig. 10 Validation of the m6A-enriched genes and candidate loci.A Validation of the m6A-enriched genes as indicated by m6A-immunoprecipitation (IP)-qPCR.B Validation of the expression levels of candidate mRNAs, lncRNAs, and circRNAs as determined by RT-qPCR

Table 1
Characteristics of the patients included in the profiling analysis

Table 2
Characteristics of the patients included in the validated cohort

Table 3
The ten most hypermethylated m6A modification peaks in CR and non-CR patients

Table 4
The ten most hypomethylated m6A modification peaks in CR and non-CR patients

Table 5
The 25 most differentially expressed mRNAs based on P value